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Abstract The statistics of the temperature anisotropies in the primordial cos- 
mic microwave background radiation field provide a wealth of information for 
cosmology and for estimating cosmological parameters. An even more acute in- 
ference should stem from the study of maps of the polarization state of the CMB 
radiation. Measuring the extremely weak CMB polarization signal requires very 
sensitive instruments. The full-sky maps of both temperature and polarization 
anisotropies of the CMB to be delivered by the upcoming Planck Surveyor satel- 
lite experiment are hence being awaited with excitement. Multiscale methods, 
such as isotropic wavelets, steerable wavelets, or curvelets, have been proposed 
in the past to analyze the CMB temperature map. In this paper, we contribute 
to enlarging the set of available transforms for polarized data on the sphere. We 
describe a set of new multiscale decompositions for polarized data on the sphere, 
including decimated and undecimated Q-U or E-B wavelet transforms and Q-U or 
E-B curvelets. The proposed transforms are invertible and so allow for applications 
in data restoration and denoising. 
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1. Introduction 



The statistical analysis of the slight intensity fluctuations in the primordial cosmic 
microwave background radiation field, for which evi dence was found for the first time in 
the early 1990's in the observations made by COBE (ISmoot et all 119921 ) . is a major issue 
in modern cosmology as these are strongly related to the cosmological scenarios describ- 
ing the properties and evolution of our Universe. In the Big Bang model, the observed 
CMB anisotropics are an imprint of primordial fiuctuations in baryon-photon density 
from a time when the temperature of the Universe was high enough above 3000 K for 
matter and radiation to be tightly coupled. At that time, the attraction of gravity and 
the repulsive radiation pressure were opposed, thus generating so-called acoustic oscil- 
lations in the baryon-photon fluid, causing peaks and troughs to appear in the power 
spectrum of the spatial anisotropics of the CMB. With the Universe cooling down as 
it expanded, matter and radiation finally decoupled. Photons were set free in a nearly 
transparent Universe, while the density fiuctuations collapsed under the effect of gravity 
into large scale structures such as galaxies or clusters of galaxies. Due to the expansion 
of the Universe, the CMB photons are now observed in the microwave range but are still 
distributed according to an almost perfect black body emission law. Another major result 
was the meas urement of the pol arization state and anisotropics of the CMB radiation 



field by DASI (|Kovac et al 



20021 ). Only a fraction of the total CMB radiation is polarized 



so that extremely sensitive instruments are needed. Polarization of the CMB radiation is 
a consequence of the Thomson scattering of photons on electrons. But for the outgoing 
population of photons to be polarized, the radiation incident on the scatterer needs to 
be anisotropic and have a quadrupolc moment. The statistics of the CMB polarization 
anisotropics are also a source of information for cosmology. Inference of cosmological pa- 
rameters from the joint statistics of the CMB anisotropics should benefit from both the 
complementarity and the redundancy of the information carried by the additional mea- 
surement of CMB polarization. Hence the full-sky maps with unprecedented sensitivity 
and angular resolution of both temperature and polarization anisotropics of the CMB 
to be delivered by the upcoming Planck Surveyor satellite experiment are awaited with 
excitement. 
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Multiscale transforms on the sphere 
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Wavelets have been us ed in many studies of CMB data ana l ysis, especially for 



Cavon et al 



Starck et al 



2001 



Holschneider 



non-Gaussianity detection ( 


Aghanim and Forni, 


I999I 


Hobson et al. 


. I999I 


Vielva et al.. 


2004; Jin et al. 


. 2005 


; Starck et al.. 


2004a 




Vielva et al. 




20061 


McEwen et al.. 


200?!. 


Continuous wavelet transforms on the sphere ( 


Antoine 




1999 


• ^ 


?enorio et al. 


1999 



19961 ) are the most used transforms. Recently, 



(120061 ) proposed an invertible isotropic undccimated wavelet transform 



(UWT) on the sphere based on spherical harmonics. A similar wavelet construction 



has been published in (jMarinucci et al 



using so-called "needlet filters". 



2008; 



Wiaux et al 



Fay and Guillouxl , 



2008 



Fay et al 



20081 ) 



20081 ) have also proposed an algorithm 



which allows us to reconstruct an image from its steerable wavelet transform. Since 
reconstruction algorithms are available, these new tools can be used for other appli- 
cations tha n non-Gaussiani t y det e ction, such as den o ising, deconvolution, com ponent 



separation (Moudden et al 



painting (lAbrial et al 



200 



2007 



1 



Bobin et al 



Abrial et al 



2008 



Delabrouille et al 



20081 ) or in- 



20081 ). Other multiscale transform s on the 



sphere such as ridgelets and curvelets have been developed (jStarck et al 



20061 ), and are 



well adapted to detect anisotropic features. It has been shown that such constructions 



are very useful for the detection of cosmic strings (jStarck et al 



Hammond et al 



2004a 



Jin et al 



2005 



20081 ) . More generally, each transform can be characterized by a given 
dictionary $ — cjix}, and the transformation consists in representing the input 

data y as a linear combination of the dictionary elements: y — — '^j^(f>kCtki where 
a are the coefficients. In the case of a wavelet transform, the dictionary is composed of 
wavelet functions, and a are the wavelet coefficients. Depending on the content of the 
data, a given dictionary may be more adapted to detect the signal of interest. A dic- 
tionary is considered as well designed for a class of signals if the transformation of any 
of these signals leads to a sparse representation, i.e. if most of the coefficients are equal 
or close to zero, while only few coefficients have a high amplitude. If the morphological 
signature of the features to be detected are not known, we should not restrict ourselves 
to analyze the data with wavelets or curvelets only, but rather use all available tools 



( Starck et al 



2004a ). 
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In this paper, we describe a set of new multiscale decompositions for polarized data 
on the sphere, including decimated and undecimated Q-U or E-B wavelet transforms, and 
Q-U or E-B curvelets. Choosing the right transform will depend on the application and 
on the typical structures occuring in the data to be analyzed. Indeed, as illustrated next, 
the different transforms are associated with different waveform dictionaries. Section [2] 
presents an orthogonal decomposition for polarized data on the sphere. Section [3] in- 
troduces a new decomposition where the modulus and the phase are processed indepen- 
dently. Sections|3]and[5]extend the isotropic wavelet transform and the curvelet transform 
on the sphere to the case of polarized data, and section [5] introduces two other wavelet 
and curvelet decompositions which are based on the E /B mode separation. The proposed 
transforms are invertible and so allow for applications in data restoration and denoising. 
Section [7] reports on denoising experiments using these new polarized transforms, thus 
demonstrating their usefulness in practice. 

2. Orthogonal representation of polarized data 

2.1. Introduction 

Figure 1. Q-U orthogonal Wavelet Transform. 

Full-sky CMB polarization data, as expected from the upcoming Planck experiment, 
consists of measurements of the Stokes parameters so that in addition to the temperature 
T map, Q and U maps are given as well. The fourth Stokes parameter commonly denoted 
y is a measure of circular polarization. In the case of CMB which is not expected to have 
circularly polarized anisotropics, V vanishes. The former three quantities, T, Q and U 
then fully describe the linear polarization state of the CMB radiation incident along 
some radial line of sight : T is the total incoming intensity, Q is the difference between 
the intensities transmitted by two perfect orthogonal polarizers the directions of which 
define a reference frame in the tangent plane, and U is the same as Q but with polarizers 
rotated 45 degrees in that tangent plane. Clearly, Q and U are not invariant through 
a rotation of angle cj) of the local reference frame around the line of sight. In fact, it is 
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easily shown that : 

Q' = cos(20)Q + sin(2(/))t/ (1) 
U' = cos{2(j))U - sm{2(j))Q 

which can also be written Q' ± ill' — e^'^'^{Q ± ill) which by definition expresses the fact 
that the quantities Q ±iU are spin-2 fields on the sphere. The suitable generalization of 
the Fourier representation for such fields is the spin-2 spherical harmonics basis denoted 
±2^m, in which we can expand : 

Q±iU = ^ ±2ae,n±2Ye,n (2) 



Figure 2. Top : examples of backprojections of Q- wavelet coefficients. Bottom : examples 
of backprojections of U-wavelet coefficients. 



2.2. Multiscale Representation 

The easiest way to build a mult i scale transform for polarized data is to use the 
Healpis 3 representation ( Gorski et al. , 2005 ) , and to apply a bi-orthogonal wavelet trans- 



form on each face of the Healpix map, separately for Q and U. Fig.[l]shows the flow-graph 
of this Q-U orthogonal wavelet transform (QU-OWT). Recall that the base resolution of 
the Healpix representation divides the sphere into twelve curvilinear quadrilateral faces of 
equal area placed on three rings around the poles and equator. Each face is subsequently 
divided into nside^ pixels of exactly equal surface but with varying shape. It follows that 
Q and U are reconstructed at position k from their wavelet coefficients w'j^, w^^, Cj ^ 
and Cjp according to : 



J 



u 



.j,pvjAp) ^j.fe(p)^ 

p j=l 



3,P 



3,P 



(3) 



^ http://healpix.jpl.nasa.go 
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which can also be written as: 



p 



J 



(4) 



p i=i 



This wavelet transform is not redundant i.e. the decomposition has the same number 
of coefficients as the input data, and it is invertible so that the Q and U maps can be 
reconstructed exactly. 

When we apply such a decomposition, we implicitly use a dictionary $ on which we 
project the data. As discussed previously, the shape of the dictionary elements, also called 
atoms, is very important to have an efficient analysis of the data. In the case of polarized 
data, it is not straightforward to imagine these shapes from Eq. [4l In order to visualize 
them, we can perform a backprojection i.e. we apply the inverse wavelet transform to 
sets of wavelet coefficients where only one coefficient is different from zero. Repeating the 
same experiment, changing only the scale and position of the non-zero coefficient allows 
us to view the different atoms in the dictionary related to the QU-OWT transform that 
we use. Fig. [5] shows examples of backprojections of Q wavelet coefficients (top) and 
backprojections of U wavelet coefficients (bottom). The shapes of the individual atoms 
do not look close to the astronomical patterns we would expect in our data. Therefore, this 
decomposition may not be optimal to analyze polarized astronomical data, although this 
would need to be confirmed in practice. The following sections describe other polarized 
wavelet transforms with different morphological properties. 

3. Module-phase non linear multiscale transform 

3.1. Introduction 

Given a polarized map in the standard Q-U representation, consider a different point 
of view and define the modulus M and phase P maps as follows : 




(5) 



Vfc, Pk = exp{i9k) where tan{9k) = Uk/Qk 



(6) 
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Because the smoothness of the Q and U maps should result in some smoothness of 
the modulus map M and the phase map P, one may consider devising a multiscale 
modulus/phase decomposition of the spin 2 field V = [Q C/]. 

The specificity of the modulus/phase decomposition of V is twofold : i) the modulus field 
is non-negative and ii ) the phase field takes its values on the unit circle 5*^. Recently, 



( Rahman et al 



2005f ) introduced a multiscale analysis technique for manifold valued 



data that will be described in the following paragraph. We then define the modulus/phase 
(MP) multiscale transform as follows : 



1. Apply a classical multiscale transform (i.e. wavelets) to the modulus map M. 

2. Apply the mult i scale analysis technique for manifold valued data described in 



I Rahman et al 



20051 ) to the phase map P. 



3.2. Decimated MP-multiscale transform 



Figures. Examples of MP-multiscale coefficients backprojection. 



Let us provide some essential notation : we assume th at the entries of the 



map P lie in a manifold A4 {e.g. A4 = S^). According to (| Rahman et al 



phase 



2005 ). take 



Pa,Pi G M and define Logpg{pi) as the log-map of pi onto the tangent space %g of Ai 
at pq. The back-projection is obtained using the inverse of the log- map Exppg. □ 
For instance, if we choose Ad = then pq — exp(i0o) and pi = exp(z0i). The Exppg 
and Logpg maps are then defined as follows : 



Vpi e S' Logpg{pi) = 6*1 -6*0 (7) 
VseM Expp^{s) = exp{i{ea + s)) (8) 

^ In differential geometry, the Exp map and Log map are generalizations of the usual expo- 
nential and logarithm function. In this paper, the manifold is a Riemannian manifold. In 
that case, the Exp map at point po, Exppg{s) is the map which takes a vector s of the tangent 
space of Ai at po and provides the point pi by travelling along the geodesic starting at po in 
the direction s. 
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The multiscale transform for manifold valued data introduced in (Rahman et al 



20051) 



is equivalent t o a two-step interpolation-refinem ent scheme similar to the lifting scheme 



described in (Daubechies and Sweldens 



19981 ). The wavelet coefficients and low pass 



approximation pixels are then computed as follows at each scale j and pixel k : 



w 



_ _ / TJf„..P 

■j+l,k 



c 



Exp,.J-Vliwf+,^,)) (10) 



The wavelet coefficient w^j^ f. at pixel k and scale j is the projection of its predic- 
tion/interpolation 'P{c^2k) onto the tangent space "^^^k+i ^ '^^2k+i- The low 
pass approximation cj^^ at scale j + 1 is computed by updating cj2k from the wavelet 
coefficient wj^^ ^. 

The main advantage of this scheme is its ability to capture local regularities while 
guaranteeing the low pass approximation to belong to the manifold M.. Indeed, the 
wavelet coefficient wf^i k pixel k and scale j + 1 is computed as the Exp map at 
S^2fe+i of an approximation :P(c|^2/c) of cj^2fe- 

Note also that even if the definitions of the Expp„ and Logpg maps involve the absolute 
phase 9{k) {i.e. tan{9{k)) — Uk/Qk), at least they only require the computation of 
differences of phases values thus avoiding the explicit manipulation of an absolute phase. 
However the non-linearity of the proposed transform is a major drawback when 
considering denoising and restoration applications. 



Illustration : 

In the case of polarized data, the entries of the phase map P lie in A4 = S^. In the 
following experiments, V and U are chosen such that : 

wf+l = ^05,P,,^,(S^2fe) (11) 
4... = ^-Pe^. (12) 

This multiscale transform is invertible and its inverse is computed as follows : 

cf^^k = Exp^.^^ ^ (^^) (13) 
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The picture in Figure [3] features some examples of backprojections of MP-multiscale 
coefficients. 



3.3. Undecimated MP-multiscale transform 



For image restoration purposes, the use of undecimated multis cale transforms has 



been shown to prov i de be tter results than decimated transforms (jStarck et al 



Starck and Murtagh 



1998; 



20061 ). The aforementioned modulus/phase multiscale analysis can 



be extended to an undecimated scheme consisting in : i) applying an undecimated wavelet 
transform to the modulus map, ii) analyzing the phase m ap P using an extensio n to the 



undecimated case of the multiscale transform described in (| Rahman et al 
case, Equations © and PH)) are replaced with the following equations : 



20051) . In that 



(15) 
(16) 



where J-{cJ ) — "YliihiLogcj ^ {cj,k-2H) with J^i^i — 1 ''^^'^ > 0- The low pass 
approximation c|^j^ ^ is then computed from a linear combination (linear filter) of a 
neighborhood {cj fc_2J/}; of weighted by the positive scalars {ft./}/. Note that from 
scale j to scale j + 1, the spatial size of the neighborhood increases by a factor 2 which 
would be equivalent to downsize by a factor 2 the band pass filter of the classical wavelet 
decomposition scheme. 



3.4. Example 

In the case of polarized data, the entries of the phase map P lie in = 5^. In the 
following experiments, J- is chosen such that : 



Exp^f, (cf+i,.) (18) 
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where : 






if 


I < 


-2 


or / > 2 


1/16 


if 


I = 


-2 


or / = 2 


1/4 


if 


I = 


-1 


or / = 1 


3/8 


if 




I 


= 



(19) 



This multiscale transform is invertible and its inverse is computed as follows : 

4,k = ^'^Pc^+,,, (-^i+i.fe) (20) 



Figure 4. Polarized field smoothing - top left : simulated synchroton emission, top right : 
same field corrupted by additive noise, bottom left : MP-multiscale reconstruction after 
setting to zero all coefficients from the three first scales, bottom right : MP-multiscale 
reconstruction after setting to zero all coefficients from the five first scales. 



Fig. |4] top shows a simulated polarized field of the synchrotron emission and its noisy 
version. We have applied the MP-multiscale transform and we remove the first three 
scales (i.e. we put all coefficients to zero) before reconstructing. The resulting image is 
shown on the bottom left of Fig. U) The bottom right of Fig. 0] corresponds to the same 
experiment, but by removing the five first scales. We can see that the field is smoother 
and smoother, but respecting the large scale structure of the field. This transform will 
be very well suited to CMB stud i es wh e re the phase is analy zed independently of the 



modulus, such as in 



Dineen et al 



( 20051) : 



Naselskv et al 



( 20051) . 



4. Polarized Wavelet Transform using Spherical Harmonics 

4.1. Isotropic Undecimated Wavelet Transform on the Sphere (UWTS) 



Figures. Q-isotropic wavelet transform backprojection (left) and U-isotropic wavelet 
backprojection (right). 
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The undccimated isotropic transform on the sphere described in (jStarck et ahl . 120061 ) 
is similar in many respects to the usual a trous isotropic wavelet transform. It is ob- 
tained using a zonal scaling function (pi^ (i?, ip) which depends only on colatitude § and 
is invariant with respect to a change in longitude Lp. It follows that the spherical har- 
monic coefficients 4>i^{l, m) of 4>i^ vanish when m ^ which makes it simple to compute 
spherical harmonic coefficients co{l,m) of cq — (j^ic * f where * stands for convolution : 



co(Z,m) = 4>i^ * f{l,m) = 



21 + 1 



(21) 



A possible scaling function (jStarck et al 



19981) . defined in the spherical harmonics rep- 



resentation, is (f>i^{l,m) — |J53(2^) where is the cubic B-spline compactly supported 
over [—2,2]. Denoting ^ rescaled version of (pi^ with cut-off frequency 2~Hc, a 

multi-resolution decomposition of / on a dyadic scale is obtained recursively : 



Co = (ki, * f 



(22) 



where the zonal low pass filters hj are defined by 



21 + 1 ' ' 



if / < and m = 



0^ (i,m) 

21 







2J+T 

otherwise 



(23) 



The cut-off frequency is reduced by a factor of 2 at each step so that in applications where 
this is useful such as compression, the numb er of samples could be reduced adequately. 
Using a pixelization scheme such as Healpix (jGorski et al.l . 120051 ) , this can easily be done 
by dividing by 2 the Healpix nside parameter when computing the inverse spherical 
harmonics transform. As in the d trous algorithm, the wavelet coefficients can be defined 
as the difference between two consecutive resolutions, Wj+i{'d, (p) — Cj{'d, tp) — Cj+i(^?, tp). 
This defines a zonal wavelet function ipi^ as 



ijjij^ {I, m) = (j) ic {I, m) — 4>ia, (/, m) 



(24) 
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Figure 6. Simulated observations on the sphere of the polarized galactic dust emission. 



With this particular choice of wavelet function, the decomposition is readily inverted 
by summing the coefficient maps on all wavelet scales 

.7 



(25) 



where we have made the simplifying assumption that / is equal to cq. Obviously, 
other wavelet functions could be used just as well, such as the needlet func- 



tion (Marinucci et al. 



2008 ) 



Figure 7. QU-Undecimated Wavelet Transform of the simulated polarized map of galac- 
tic dust emission shown in figure 



4.2. Extension to Polarized Data 

By applying the above scalar isotropic wavelet transform to each component T, Q, U 
of a polarized map on the sphere, we have : 

T{§,^) =cJ(i?,^)H-E/=i«^J(^,^) (26) 
Q(^,^) =c;3(^,^)-hE/=i<(^,^) 

where stands for the low resolution approximation to component X and Wj^ is the 
map of wavelet coefficients of that component on scale j. This leads to the following 
decomposition : 

(Q ± ^Um = (c? ± c%)[k] + J2iwf ± 0[fc] (27) 

7 = 1 

Fig [5] shows the backprojection of a Q-wavelet coefficient (left) and a [/-wavelet coef- 
ficient (right). Fig. [7] shows the undecimated isotropic polarized wavelet transform of 
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the dust image shown on Fig. [5] using six scales, i.e. five wavelet scales and the coarse 
approximation. 



5. Polarized Curvelet Transform 



Figures. Top, Q-curvelet backprojection (left) and zoom (right). Bottom, U-curvelet 
backprojection (left) and zoom. 



The 2D ridgelet transform (jCandes and Donohd . Il999ar ) was developed in an at 



tempt to overcome some limitations inherent in former multiscale methods e.g. the 
2D wavelet, when handling smooth images with edges i.e. singularities along smooth 
curves. Ridgelets are translation invariant ridge functions with a wavelet profile in 
the normal direction. Although ridgelets provide sparse representations of smooth im- 
ages with straight edges, they fail to efiiciently handle edges along curved lines. This 
is the framework for curvel ets which were given a first mathematical description in 



( Candes and Donoho 



IQQQbT ). Basically, the curvelet dictionary is a multiscale pyramid 



of localized directional functions with anisotropic support obeying a specific parabolic 
scaling such that at scale 2~^, its length is 2^^/^ and its width is This is motivated 
by the parabolic scaling property of smooth curves. Other properties of the curvelet 
tra nsform as well as decisive opt imality results in approximation theory are reported 



in ( Candes and Donoho 



sive op 
1999bic 



dj). Notably, curvelets provide optimally sparse repre- 



sentations of manifolds which are smooth a way from edge singul a rities along smooth 



curves. Several dig i tal cu rvelet transforms (jDonoho and Duncan 



2002 : 



Candes et al 



200C; 



Starck et al 



20061 ) have been proposed which attempt to preserve the essential 



properties of the continuous curvelet trans form and many papers (jStarck et al 



Herrmann et al 



2008 



Starck et al 



2004b : 



2004br ) report on their successful application in im- 



age processing experiments. The so-called first generation discrete curvelet described in 



( Donoho and Duncan 



2000 



Starck et al 



20021 ) consists in applying the ridgelet trans- 



form to sub- images of a wavelet decomposition of the original image. By construction, the 
sub-images are well localized in space and frequency and the subsequent ridgelet trans- 
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form provides the necessary directional sensitivity. This latter implementation in combi- 
nation with the good geometric properties of the Healpix pixeliz ation scheme, inspired 
the digital curvelet transform on the sphere (|Starck et al.l . 120061 ). The digital curvelet 
transform on the sphere is clearly invertible in the sense that each step of the overall 
transform is itself invertible. The curvelet transform on the sphere has a redundancy 
factor of 16 J + 1 when J scales are used, which may be a problem for handling huge 
data sets such as from the future Planck-Surveyor experiment. This can be reduced by 
substituting the pyramidal wavelet transform to the undecimated wavelet transform in 
the above algorithm. More details on the wavelet, ridgelet, curvelet algorithms on the 



sphere can be found in (jStarck et al 



20061 ) and software related to these new transforms 



is available from the web page http://jstarck.free.fr/mrs.htmli As for the isotropic wavelet 
on the sphere, a straightforward extension to polarized data will consist in applying suc- 
cessively the curvelet transform on the sphere to the three components T, Q and U. 
Figure [8] shows the backprojection of a Q-curvelet coefficient and U-curvelet coefficient. 
Clearly, the shapes of these polarized curvelet functions are very different from the polar- 
ized wavelet functions. In the next section, we will build very different dictionaries using 
the E/B mode decomposition. 



6. Polarized E/B Wavelet and E/B Curvelet 

6.1. Introduction 

We have seen that the generalization of the Fourier representation for polarized data 
on the sphere is the spin-2 spherical harmonics basis denoted ±2Y£rn' 



Q ±iU = ^ ±2a£m±2Yl 



im 



(28) 



At this point, it is convenient (jZaldarriagal 119981 ) to introduce the two quantities 
denoted E and B which are defined on the sphere by 



lYem — Ef m ~^(2a^m + ~2aem)Yi, 



(29) 



where stands for the usual spin spherical harmonics basis functions. The quantities 
E and B are derived by applying the spin lowering operator twice to Q + iU and the spin 
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raising operator twice to Q — ill so that E and B are real scalar fields on the sphere, 
invariant through rotations of the local reference frame. The normalization of af^ and 
^fm cho sen in the latter defini t ion is purely conventiona l but it appears to be rather 



Zaldarriaga and SeljakI (119971 ): 



Bunnetal 



([20031). StiU, we could multiply a 



popular 

and af^ by some Af^ and we would have just as good a representation of the initial po- 
larization maps. Through a change of parity E will remain invariant whereas the sign 
of pseudo-scalar B will change. The E and B modes defined here are not so different 
from the gradient (i.e. curl free) and curl {i.e. divergence free) components encountered 
in the analysis of vector fields. Finally, the spatial anisotropics of the Gaussian CMB 
temperature and polarization fields are completely characterized in this new linear rep- 
resentation by the power spectra and cross spectra of T, E and B. Thanks to the different 
parities of T and E on one side and B on the other, the sufficient statistics reduce to 
only four spectra namely Cf^,Cf^,Cf^,Cf^. For a given cosmological model, it is 
possible to give a theoretical prediction of these spectra. Aiming at inverting the model 
and inferring the cosmological parameters, an important goal of CMB temperature and 
polarization data analysis is then to estimate the latter power spectra, based on sampled, 
noisy sometimes incomplete T, Q, [/ spherical maps. 



6.2. E/B Isotropic Wavelet 



Following the above idea of representing CMB polarization maps by means of E and 
B modes, we propose a formal extension of the previous undecimated isotropic wavelet 
transform that will allow us to handle linear polarization data maps T,Q,U on the 
sphere. Practically, the maps we consider are pixelized using for instance the Healpix 
pixelization scheme. In fact, we are not concerned at this point with the recovery of E 
and B modes from pixelized or incomplete data maps which itself is not a trivial task. 
The extension of the wavelet transform on the sphere we describe here makes use of 
the E and B representation of polarized maps described above in a formal way. Given 
polarization data maps T, Q and [/, the proposed wavelet transform algorithm consists 
of the following steps : 
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Figure 9. top : Q and U CMB polarization data maps from channel ka of the WMAP 
experiment, left : low pass and wavelet coefficients in three scales of the formal E mode, 
right : low pass and wavelet coefficients in three scales of the formal B mode. 



1. Apply the spin ±2 spherical harmonics transform to Q + iU and Q — ill. 
Practically, the Healpix software package provides an implementation of this 
transform for maps that use this pix elization scheme. Ot herwise, a fast imple- 
mentation was recently proposed by 
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2. Combine the decomposition coefficients 2<iem and _2a^m from the first step 
into af^ and a£„ and build formal E and B maps associated with Q and U by 
applying the usual inverse spherical harmonics transform, as in equation 1291 
For numerical and algorithmic purposes, it may be efficient to stay with the 
spherical harmonics representation of E and B. 

3. Apply the undecimated isotropic transform on the sphere described above to 
map T and to the E, B representation of the polarization maps. 

The wavelet coefficient maps wj, ^ and the low resolution approximation maps cj, 
Cj , Cj are obtained by applying the isotropic undecimated wavelet transform described 
in section 14. f I to the T, E, B representation of the polarized data. Figure [9 shows the 
result of applying the proposed transform to the polarized CMB data map ka f| from the 
WMAP experiment. The top two images show the initial Q and U maps while the subse- 
quent maps are the low pass and wavelet coefficients' maps in a four scale decomposition. 
The scaling function we used is a cubic box spline as proposed in section [4T] The wavelet 
coefficients were obtained as the difference between two successive low pass approxima- 
tions of the multiresolution decomposition of the E and B maps. The proper choice for 
the scaling and wavelet functions will depend on the application and the existence of 
constraints to be enforced. 



available at http:/ /lambda. gsfc.nasa.gov/product /map /current / 1 



Reconstruction 
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Figure 10. E-isotropic wavelet transform backprojection (left) and B-isotropic wavelet 
backprojection (right). 

Obviously, the transform described above is invertible and the inverse transform is 
readily obtained by applying the inverse of each of the three steps in reverse order. If, as 
in the example decomposition above, we take the wavelet function to be the difference 
between two successive low pass approximations, the third step is inverted by simply 
summing the last low pass approximation with the maps of wavelet coefficients from all 
scales as in equation [55] : 



J 




J 



E 



J 




(30) 



where Cj" stands for the low resolution approximation to component X and w-^ is the 
map of wavelet coefficients of that component on scale j. Finally, noting that : 




ra 





(31) 



1 





U 



rn 




(32) 
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the initial representation of the polarized data in terms of T, Q and U maps is recon- 
structed from its wavelet coefficients using the following equations : 

(33) 



T = 






Q = cf+ 


+ IC j ■ 


- , v^J f E. 


u = cr 


■ E,- 
-ICj 


- , v^J f B, 



■ B. 



■ E,- 



] 

where 

-cf^yi„Z+„ and cJ-=cJ5]r/„Z,-^ (34) 



with W'^ denoting the transpose conjugate of W so that WW'' is the scalar dot product 
of W and W while W^W is an operator (or matrix) acting on its left hand side as a 
projection along W and reconstruction along W. In practice, the Healpix software package 
provides us with an implementation of the forward and inverse spin and spin 2 spherical 
harmonics transforms which we need to implement the proposed inverse transform given 
by equations [33] and [34l Clearly, as mentioned earlier in section llTl] we could have chosen 
some other wavelet function than merely the difference between two consecutive scaling 
functions, and the transformation would still be nearly as simple to invert. Fig 1 101 shows, 
on the left, backprojections of E-wavelet coefficients, and, on the right, backprojections 
of B-wavelet coefficients on the right hand side at different scales. 

E-B Curvelet 



Figure 11. E-curvelet coefficient backprojection. 



Figure 12. B-curvelet coefficient backprojection. 

Similarly to the EB-wavelet constructions, we can easily construct an EB-curvelet 
transform by first computing the E and B components using the spin ±2 spherical har- 
monics transform, and then applying a curvelet transform on the sphere separately on 
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Table 1. Error (in dB) for both the Q and U components between the original dust 
image and respectively its noisy version and the results obtained by hard thresholding 
representations of the noisy data in different dictionaries. 



Q Component 




Noisy image 


20.01 


13.98 


6.03 


-0.007 


-6.02 


QU-UWT 


34.91 


32.41 


28.88 


26.41 


23.68 


EB-UWT 


33.52 


30.63 


27.83 


25.86 


23.91 


QU-CUR 


34.14 


31.12 


28.56 


26.58 


24.56 


EB-CUR 


33.36 


30.60 


28.17 


26.30 


23.99 


OWT 


34.07 


30.77 


27.08 


24.57 


21.33 


Mod-phase 


30.51 


26.42 


20.97 


16.25 


11.05 


U Component 




Noisy image 


19.99 


13.97 


6.01 


-0.001 


-6.01 


QU-UWT 


40.88 


39.22 


36.77 


35.50 


31.90 


EB-UWT 


38.80 


36.71 


35.41 


34.67 


32.23 


QU-CUR 


39.54 


37.85 


36.68 


35.83 


32.38 


EB-CUR 


39.64 


37.72 


35.33 


33.51 


29.30 


OWT 


39.18 


37.29 


33.30 


29.06 


23.56 


Mod-phase 


30.76 


26.81 


21.31 


16.83 


11.36 



each of these two components. Fig. [TT] shows the backprojection of an E-curvelet coeffi- 
cient and Fig. [T2l shows the backprojectionof a B-curvelet coefficient. 



Figure 13. Top, simulated input polarized image (left) and noisy polarized field (right) 
Bottom, denoising of the polarized field using the EB-isotropic undecimated wavelets 
(left) and the EB-curvelets (right). 



20 Starck et ai: Polarized wavelets and curvelets 

Table 2. Error (in dB) for both the Q and U components between the original syn- 
chrotron image and respectively its noisy version and the results obtained by the hard 
thresholding using different dictionaries. 



Q Component 




Noisy image 


19.99 


13.98 


6.02 


-0.005 


-6.02 


QU-UWT 


33.40 


31.27 


27.78 


24.71 


21.16 


EB-UWT 


33.03 


30.53 


26.79 


23.86 


20.76 


QU-CUR 


35.04 


32.92 


29.19 


26.19 


22.37 


EB-CUR 


34.83 


32.09 


27.89 


24.81 


21.06 


OWT 


33.20 


30.10 


25.72 


22.62 


19.39 


MP-UWT 


26.76 


23.05 


17.85 


13.69 


8.97 


U Component 




Noisy image 


19.99 


13.97 


6.01 


-0.007 


-6.01 


QU-UWT 


32.90 


31.25 


29.04 


27.51 


25.27 


EB-UWT 


33.22 


31.43 


29.11 


26.99 


24.35 


QU-CUR 


33.79 


31.75 


29.28 


27.71 


25.46 


EB-CUR 


34.05 


31.98 


29.38 


27.17 


24.07 


OWT 


32.69 


30.49 


27.87 


25.49 


21.87 


MP-UWT 


26.76 


22.76 


17.85 


13.86 


9.41 



7. Experiments : Application to denoising 

Thanks to the invertibility of the different proposed transforms for polarization maps 
on the sphere, it is possible to use these transformations for restoration applications. The 
denoising algorithm we use here consists in three consecutive steps : 

1. Apply the chosen polarized multiscale transform. 

2. Set to zero those coefficients whose absolute values are below a given threshold. We 
have used a threshold equal to five times the noise standard deviation. 

3. Reconstruct the denoised field using the inverse transform. 



More sophisticated thresholding strategies exist (jStarck and Murtaghl . 120061) which can 
be used just as well on coefficients of polarized wavelets and curvelets. 



Starck et at: Polarized wavelets and curvelets 21 

Figure [Ol illustrates the results of a simple denoising experiment. Noise was added 
to a simulated dust image (see Fig. [T51 top left and right), and the noisy QU-field was 
filtered by thresholding either the EB-isotropic wavelet coefhcients of the polarized dust 
map (Fig. [13] bottom left) or the EB-isotropic curvelet coefficients (Fig. [13] bottom right). 
Both decompositions produce nice visual results. 

In order to be more quantitative, we used two test images (synchrotron and dust) with 
different noise levels. The noise levels were taken equal to 0.1, 0.2, 0.5, 1. and 2 times the 
standard deviation of the noise-free image. On each noisy image, we applied six different 
transformations, thresholded the coefficients and reconstructed the filtered images. The 
transforms we used are the QU- wavelets (QU-UWT), the EB-wavelets (EB-UWT), the 
QU-curvelet (QU-CUR), the EB-curvelet (EB-CUR), the biorthogonal wavelet transform 
(OWT) and the modulus-phase undecimated multiscale transform (MP-UWT). For each 
filtered image Q (resp. J7), we computed the Signal-to-Noise Ratio (SNR) in dB : 

E=lQ\og,^{a]/cyl) (35) 

where cr/ is the standard deviation of the noise free original image component Q (resp. 
U) and Ue is the standard deviation of the error image i.e. the difference between the 
noise-free component Q (resp. U) and the filtered component Q (resp. U). These errors 
are given in Table ([T]) for the dust image and in Table ^ for the synchrotron image. 

It is clear from the above results that the different transforms described here do not 
perform as well on this specific numerical experiment. For the dust, QU- wavelets are 
better when the noise level is not so high, while curvelets do better when the noise and 
signal levels are of the same order. This can be explained by the fact that structures on 
small scales are more or less isotropic, and therefore better represented by wavelets, while 
large scale structures are anisotropic and therefore better analyzed using curvelets. When 
increasing the noise level, structures on the smallest scales are no longer detected by either 
of the two dictionaries. Only large scale features are detectable, and curvelets do this job 
better. Dealing with the polarized synchrotron map, curvelets do better than wavelets in 
all cases experimented with here. Although the bi-orthogonal wavelet transform is clearly 
not as good as the others in these experiments, it could nevertheless be very useful in 
situations where computation time is an important issue. Indeed, since it doesn't make 
use of the spherical harmonics transform and also because it is not redundant, it is a 
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very fast transform. Finally the worse results were obtained using the Modulus-Phase 
multiscale transform. This can be explained by the fact that the Gaussian approximation 
we made for the noise in the wavelet transform of the modulus is not accurate enough. 
Furthermore, it is not clear what is the best way to correct the phase wavelet coefficients 
from the noise. A better understanding of the noise behavior after transformation is 
clearly required before the Mod-phase multiscale transform can be used for restoration 
purposes. However, for other applications such as non-Gaussianity studies, the latter 
transform could prove an interesting tool to use as well. 

8. Conclusion 

The present contribution has enlarged the set of available representations of polar- 
ized data on the sphere. We described the construction of new multiscale decompositions, 
namely the modulus-phase transform, the QU and EB isotropic undecimated wavelet and 
curvelet transforms for polarized data on the sphere. The latter are derived as extensions 
of the undecimated wavelet and curvelet transforms for scalar maps on the sphere. The 
proposed extensions use a formal representation of T, Q and U polarization data maps 
in terms of the scalar T, E and pseudo-scalar B maps. The proposed transforms are 
invertible allowing for applications in image restoration and denoising as was shown in a 
preliminary experiment. Ongoing research is concerned with assessing and understanding 
the merits of the different transforms we described, in such problems as restoring, de- 
noising or inpainting of sparse linearly polarized data, as well as in the blind separation 
of mixed linearly polarized components. 

Software 

The software related to this paper, MRS/POL, and its full documentation will be 
included in the next version of the the MRS software which is available from the following 
web page: 

\protect\vrule widthOpt\protect\href {http : //j starck . free . f r/mrs .html}{http : //j starck. free . f r/mrs .html} 



Acknowledgements. Some of the results in this paper have been derived using the Healpix 
i Gorski et aLlboOsI ). 
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